We study mainly algorithms for real symmetric matrices, which are most commonly used in the applications described in this course.
For more details, see I. Slapničar, Symmetric Matrix Eigenvalue Techniques and the references therein.
The reader should be familiar with basic linear algebra concepts and facts on eigenvalue decomposition and perturbation theory
The reader should be able to apply adequate algorithm to a given problem, and to assess accuracy of the solution.
If the value of a function $f(x)$ is computed with an algorithm $\mathrm{alg(x)}$, the algorithm error is
$$ \|\mathrm{alg(x)}-f(x)\|, $$and the relative algorithm error is
$$ \frac{\| \mathrm{alg}(x)-f(x)\|}{\| f(x) \|}, $$in respective norms. Therse errors can be hard or even impossible to estimate directly.
In this case, assume that $f(x)$ computed by $\mathrm{alg}(x)$ is equal to exact value of the function for a perturbed argument,
$$ \mathrm{alg}(x)=f(x+\delta x), $$for some backward error $\delta x$.
Algoritam is stable is the above equality always holds for small $\delta x$.
The eigenvalue decomposition (EVD) of a real symmetric matrix $A=[a_{ij}]$ is $A=U\Lambda U^T$, where $U$ is a $n\times n$ real orthonormal matrix, $U^TU=UU^T=I_n$, and $\Lambda=\mathop{\mathrm{diag}} (\lambda_1,\dots, \lambda_n)$ is a real diagonal matrix.
The numbers $\lambda_i$ are the eigenvalues of $A$, the vectors $U_{:i}$, $i=1,\dots,n$, are the eigenvectors of $A$, and $AU_{:i}=\lambda_i U_{:i}$, $i=1,\dots,n$.
If $|\lambda_1|> |\lambda_2| \geq \cdots \geq |\lambda_n|$, we say that $\lambda_1$ is the dominant eigenvalue.
Deflation is a process of reducing the size of the matrix whose EVD is to be determined, given that one eigenvector is known.
The shifted matrix of the matrix $A$ is the matrix $A-\mu I$, where $\mu$ is the shift.
Power method starts from unit vector $x_0$ and computes the sequences $$ \nu_k=x_k^T A x_k, \qquad x_{k+1}= A x_k / \| A x_k \|, \qquad k=0,1,2,\dots, $$ until convergence. Normalization of $x_k$ can be performed in any norm and serves the numerical stability of the algorithm (avoiding overflow or underflow).
Inverse iteration is the power method applied to the inverse of a shifted matrix: $$ \nu_k=x_k^T A x_k, \quad x_{k+1}= (A-\mu I)^{-1} x_k, \quad x_{k+1} = x_{k+1}/\|x_{k+1}\|, \quad k=0,1,2,\dots. $$
QR iteration starts from the matrix $A_0=A$ and forms the sequence of matrices $$ A_k=Q_kR_k \quad \textrm{(QR factorization)}, \qquad A_{k+1}=R_kQ_k,\qquad k=0,1,2,\ldots $$
Shifted QR iteration is the QR iteration applied to a shifted matrix: $$ A_k-\mu I=Q_kR_k \quad \textrm{(QR factorization)}, \quad A_{k+1}=R_kQ_k+\mu I ,\quad k=0,1,2,\ldots $$
If $\lambda_1$ is the dominant eigenvalue, $x_0$ is not orthogonal to $U_{:1}$, and the norm is $2$-norm, then $\nu_k\to \lambda_1$ and $x_k\to U_{:1}$. In other words, the power method converges to the dominant eigenvalue and its eigenvector.
The convergence is linear in the sense that $$ |\lambda_1-\nu_k|\approx \left|\frac{c_2}{c_1}\right| \left| \frac{\lambda_2}{\lambda_1}\right|^k,\qquad \|U_{:1}-x_k\|_2 =O\bigg(\bigg| \frac{\lambda_2}{\lambda_1}\bigg|^k\bigg)\!, $$ where $c_i$ is the coefficient of the $i$-th eigenvector in the linear combination expressing the starting vector $x_0$.
Since $\lambda_1$ is not available, the convergence is determined using residuals: if $\| Ax_k-\nu_k x_k\|_2\leq tol$, where $tol$ is a user prescribed stopping criterion, then $|\lambda_1-\nu_k|\leq tol$.
After computing the dominant eigenpair, we can perform deflation to reduce the given EVD for $A$ to the one of size $n-1$ for $A_1$: $$ \begin{bmatrix} U_{:1} & X \end{bmatrix}^T A \begin{bmatrix} U_{:1} & X \end{bmatrix}= \begin{bmatrix} \lambda_1 & \\ & A_1 \end{bmatrix}, \quad \begin{bmatrix} U_{:1} & X \end{bmatrix} \textrm{orthonormal}, \quad A_1=X^TAX. $$
The EVD of the shifted matrix $A-\mu I$ is $U(\Lambda-\mu I) U^T$.
Inverse iteration requires solving the system of linear equations $(A-\mu I)x_{k+1}= x_k$ for $x_{k+1}$ in each step. At the beginning, LU factorization of $A-\mu I$ needs to be computed, which requires $2n^3/3$ operations. In each subsequent step, two triangular systems need to be solved, which requires $2n^2$ operations.
If $\mu$ is close to some eigenvalue of $A$, the eigenvalues of the shifted matrix satisfy $|\lambda_1|\gg |\lambda_2|\geq\cdots\geq |\lambda_n|$, so the convergence of the inverse iteration method is fast.
If $\mu$ is very close to some eigenvalue of $A$, then the matrix $A-\mu I$ is nearly singular, so the solutions of linear systems may have large errors. However, these errors are almost entirely in the direction of the dominant eigenvector so the inverse iteration method is both fast and accurate!!
We can further increase the speed of convergence of inverse iterations by substituting the shift $\mu$ with the Rayleigh quotient $\nu_k$ at the cost of computing new LU factorization.
Matrices $A_k$ and $A_{k+1}$ from both QR iterations are orthogonally similar, $A_{k+1}=Q_k^TA_kQ_k$.
The QR iteration method is essentially equivalent to the power method and the shifted QR iteration method is essentially equivalent to the inverse power method on the shifted matrix.
The straightforward application of the QR iteration requires $O(n^3)$ operations per step, so better implementation is needed.
In order to keep the programs simple, in the examples below we do not compute full matrix of eigenvectors.
In [1]:
function myPower(A::Matrix,x::Vector,tol::Real)
y=A*x
ν=x⋅y
steps=1
while norm(y-ν*x)>tol
x=y/norm(y)
y=A*x
ν=x⋅y
# display(x)
steps+=1
end
ν, y/norm(y), steps
end
Out[1]:
In [2]:
import Random
Random.seed!(421)
using LinearAlgebra
n=6
A=Matrix(Symmetric(rand(-9:9,n,n)))
Out[2]:
In [3]:
ν,x,steps=myPower(A,ones(n),1e-10)
Out[3]:
In [4]:
λ=eigvals(A)
Out[4]:
In [5]:
# Speed of convergence
(18/21)^100
Out[5]:
In [7]:
findmax(broadcast(abs,λ))
Out[7]:
In [8]:
k=findmax(broadcast(abs,λ))[2]
Out[8]:
In [9]:
[eigvecs(A)[:,k] x]
Out[9]:
In [10]:
# Deflation
function myDeflation(A::Matrix,x::Vector)
X,R=qr(x)
# To make sure the returned matrix symmetric use
# full(Symmetric(X[:,2:end]'*A*X[:,2:end]))
# display(X'*A*X)
X[:,2:end]'*A*X[:,2:end]
end
Out[10]:
In [11]:
A₁=myDeflation(A,x)
Out[11]:
In [12]:
eigvals(A₁)
Out[12]:
In [13]:
myPower(A₁,ones(size(A₁,1)),1e-10)
Out[13]:
In [14]:
# Put it all together - eigenvectors are omitted for the sake of simplicty
function myPowerMethod(A::Matrix{T}, tol::Real) where T
n=size(A,1)
S = T==Int ? Float64 : T
λ=Vector{S}(undef,n)
for i=1:n
λ[i],x,steps=myPower(A,ones(n-i+1),tol)
A=myDeflation(A,x)
end
λ
end
Out[14]:
In [15]:
myPowerMethod(A,1e-10)
Out[15]:
In [16]:
# QR iteration
function myQRIteration(A::Matrix, tol::Real)
steps=1
while norm(tril(A,-1))>tol
Q,R=qr(A)
A=R*Q
steps+=1
end
A,steps
end
Out[16]:
In [17]:
B,steps=myQRIteration(A,1e-5)
B
Out[17]:
In [18]:
steps
Out[18]:
In [19]:
diag(B)
Out[19]:
The following implementation of $QR$ iteration requires a total of $O(n^3)$ operations:
One step of QR iterations on $T$ requires $O(n)$ operations if only $\Lambda$ is computed, and $O(n^2)$ operations if $Q$ is accumulated, as well.
Given vector $v$, Householder reflector is a symmetric orthogonal matrix
$$ H= I - 2\, \frac{v v^T}{v^Tv}. $$Given $c=\cos\varphi$, $s=\sin\varphi$, and indices $i<j$, Givens rotation matrix is an orthogonal matrix $G(c,s,i,j)$ which is equal to the identity matrix except for elements
$$ G_{ii}=G_{jj}=c,\qquad G_{ij}=-G_{ji}=s. $$Given vector $x$, choosing $v=x+\mathrm{\mathop{sign}}(x_{1})\,\|x\|_2\, e_1$ yields the Householder reflector which performs the QR factorization of $x$: $$ Hx=-\mathop{\mathrm{sign}}(x_{1})\, \|x\|_2\, e_1. $$
Given a 2-dimensional vector $\begin{bmatrix}x\\y\end{bmatrix}$, choosing $r= \sqrt{x^2+y^2}$, $c=\frac{x}{r}$ and
$s=\frac{y}{r}$, gives the Givens roatation matrix such that
$$
G(c,s,1,2)\cdot \begin{bmatrix}x\\y\end{bmatrix}=\begin{bmatrix}r\\ 0 \end{bmatrix}.
$$
The hypotenuse $r$ is computed using the hypot() function in order to avoid underflow or overflow.
Tridiagonal form is not unique.
The reduction of $A$ to tridiagonal matrix by Householder reflections is performed as follows. Let $$ A=\begin{bmatrix} \alpha & a^T \\ a & B \end{bmatrix}. $$ Let $v=a+\mathrm{\mathop{sign}}(a_{1})\,\|a\|_2\, e_1$, let $H$ be the corresponding Householder reflector and set $$ H_1=\begin{bmatrix} 1 & \\ & H \end{bmatrix}. $$ Then $$ H_1AH_1= \begin{bmatrix} \alpha & a^T H \\ Ha & HBH \end{bmatrix} = \begin{bmatrix} \alpha & \nu e_1^T \\ \nu e_1 & A_1 \end{bmatrix}, \quad \nu=-\mathop{\mathrm{sign}}(a_{1})\, \|a\|_2. $$ This step annihilates all elements in the first column below the first subdiagonal and all elements in the first row to the right of the first subdiagonal. Applying this procedure recursively yields the tridiagonal matrix $T=X^T AX$, where $X=H_1H_2\cdots H_{n-2}$.
$H$ does not depend on the normalization of $v$. With the normalization $v_1=1$, $a_{2:n-1}$ can be overwritten by $v_{2:n-1}$, so $v_1$ does not need to be stored.
$H$ is not formed explicitly - given $v$, $B$ is overwritten with $HBH$ in $O(n^2)$ operations by using one matrix-vector multiplication and two rank-one updates.
When symmetry is exploited in performing rank-2 update, tridiagonalization requires $4n^3/3$ operations. Instead of performing rank-2 update on $B$, one can accumulate $p$ transformations and perform rank-$2p$ update. This block algorithm is rich in matrix--matrix multiplications (roughly one half of the operations is performed using BLAS 3 routines), but it requires extra workspace for $U$ and $V$.
If $X$ is needed explicitly, it can be computed from the stored Householder vectors $v$. In order to minimize the operation count, the computation starts from the smallest matrix and the size is gradually increased: $$ H_{n-2},\quad H_{n-3}H_{n-2},\dots,\quad X=H_1\cdots H_{n-2}. $$ A column-oriented version is possible as well, and the operation count in both cases is $4n^3/3$. If the Householder reflectors $H_i$ are accumulated in the order in which they are generated, the operation count is $2n^3$.
The backward error bounds for functions myTridiag() and myTridiagX() are as follows:
The computed matrix $\tilde T$ is equal to the matrix which
would be obtained by exact tridiagonalization of some perturbed matrix $A+\Delta A$,
where $\|\Delta A\|_2 \leq \psi \varepsilon \|A\|_2$ and $\psi$ is a
slowly increasing function of $n$.
The computed matrix $\tilde X$ satisfies $\tilde X=X+\Delta X$, where
$\|\Delta X \|_2\leq \phi \varepsilon$
and $\phi$ is a slowly increasing function of $n$.
Tridiagonalization using Givens rotations requires $\frac{(n-1)(n-2)}{2}$ plane rotations, which amounts to $4n^3$ operations if symmetry is properly exploited. The operation count is reduced to $8n^3/3$ if fast rotations are used. Fast rotations are obtained by factoring out absolutely larger of $c$ and $s$ from $G$.
Givens rotations in the function myTridiagG() can be performed in different
orderings. For example, the elements in the first column and row can be annihilated by
rotations in the planes $(n-1,n)$, $(n-2,n-1)$, $\ldots$, $(2,3)$.
Givens rotations act more selectively than Householder reflectors, and
are useful if $A$ has some special structure, for example, if $A$ is a banded matrix.
Error bounds for function myTridiagG() are the same as above,
but with slightly different functions $\psi$ and $\phi$.
The block version of tridiagonal reduction is implemented in the LAPACK subroutine DSYTRD. The computation of $X$ is implemented in the subroutine DORGTR. The size of the required extra workspace (in elements) is $lwork=nb*n$, where $nb$ is the optimal block size (here, $nb=64)$, and it is determined automatically by the subroutines. The subroutine DSBTRD tridiagonalizes a symmetric band matrix by using Givens rotations. There are no Julia wappers for these routines yet!
In [20]:
function myTridiag(A::Matrix)
# Normalized Householder vectors are stored in the lower
# triangular part of A below the first subdiagonal
n=size(A,1)
T=Float64
A=map(T,A)
v=Vector{T}(undef,n)
Trid=SymTridiagonal(zeros(n),zeros(n-1))
for j = 1 : n-2
μ = sign(A[j+1,j])*norm(A[j+1:n, j])
if μ != zero(T)
β =A[j+1,j]+μ
v[j+2:n] = A[j+2:n,j] / β
end
A[j+1,j]=-μ
A[j,j+1]=-μ
v[j+1] = one(T)
γ = -2 / (v[j+1:n]⋅v[j+1:n])
w = γ* A[j+1:n, j+1:n]*v[j+1:n]
q = w + γ * v[j+1:n]*(v[j+1:n]⋅w) / 2
A[j+1:n, j+1:n] = A[j+1:n,j+1:n] + v[j+1:n]*q' + q*v[j+1:n]'
A[j+2:n, j] = v[j+2:n]
end
SymTridiagonal(diag(A),diag(A,1)), tril(A,-2)
end
Out[20]:
In [21]:
A
Out[21]:
In [22]:
T,H=myTridiag(A)
Out[22]:
In [23]:
T
Out[23]:
In [24]:
[eigvals(A) eigvals(T)]
Out[24]:
In [25]:
H
Out[25]:
In [26]:
v1=[1;H[3:6,1]]
Out[26]:
In [27]:
H0=I-2*v1*v1'/(v1'*v1)
Out[27]:
In [28]:
H1=[1 zeros(1,5);zeros(5,1) H0]
Out[28]:
In [29]:
H1*A*H1
Out[29]:
In [30]:
# Extract X
function myTridiagX(H::Matrix)
n=size(H,1)
T=Float64
X = Matrix{T}(I,n,n)
v=Vector{T}(undef,n)
for j = n-2 : -1 : 1
v[j+1] = one(T)
v[j+2:n] = H[j+2:n, j]
γ = -2 / (v[j+1:n]⋅v[j+1:n])
w = γ * X[j+1:n, j+1:n]'*v[j+1:n]
X[j+1:n, j+1:n] = X[j+1:n, j+1:n] + v[j+1:n]*w'
end
X
end
Out[30]:
In [31]:
X=myTridiagX(H)
Out[31]:
In [32]:
# Fact 7: norm(ΔX)<ϕ*eps()
X'*X
Out[32]:
In [33]:
X'*A*X
Out[33]:
In [34]:
# Tridiagonalization using Givens rotations
function myTridiagG(A::Matrix)
n=size(A,1)
T=Float64
X=Matrix{T}(I,n,n)
for j = 1 : n-2
for i = j+2 : n
G,r=givens(A,j+1,i,j)
A=(G*A)*G'
# display(A)
X*=G'
end
end
SymTridiagonal(diag(A),diag(A,1)), X
end
Out[34]:
In [35]:
# Let us take a look at the `givens()` functions
methods(givens)
Out[35]:
In [36]:
?givens
Out[36]:
In [37]:
Tg,Xg=myTridiagG(map(Float64,A))
Out[37]:
In [38]:
Tg
Out[38]:
In [39]:
Xg'*Xg
Out[39]:
In [40]:
Xg'*A*Xg
Out[40]:
In [41]:
T
Out[41]:
In [42]:
Q,R=qr(Matrix(T))
Out[42]:
In [43]:
R*Q
Out[43]:
Let $T$ be a real symmetric tridiagonal matrix of order $n$ and $T=Q\Lambda Q^T$ be its EVD.
Each step of the shifted QR iterations can be elegantly implemented without explicitly computing the shifted matrix $T-\mu I$.
Wilkinson's shift $\mu$ is the eigenvalue of the bottom right $2\times 2$ submatrix of $T$, which is closer to $T_{n,n}$.
The stable formula for the Wilkinson's shift is $$ \mu=T_{n,n}- \frac{T_{n,n-1}^2}{\tau+\mathop{\mathrm{sign}}(\tau)\sqrt{\tau^2+T_{n,n-1}^2}},\qquad \tau=\frac{T_{n-1,n-1}-T_{n,n}}{2}. $$
Wilkinson's shift is the most commonly used shift. With Wilkinson's shift, the algorithm always converges in the sense that $T_{n-1,n}\to 0$. The convergence is quadratic, that is, $|[T_{k+1}]_{n-1,n}|\leq c |[T_{k}]_{n-1,n}|^2$ for some constant $c$, where $T_k$ is the matrix after the $k$-th sweep. Even more, the convergence is usually cubic. However, it can also happen that some $T_{i,i+i}$, $i\neq n-1$, becomes sufficiently small before $T_{n-1,n}$, so the practical program has to check for deflation at each step.
Chasing the Bulge. The plane rotation parameters at the start of the sweep are computed as if the shifted $T-\mu I$ has been formed. Since the rotation is applied to the original $T$ and not to $T-\mu I$, this creates new nonzero elements at the positions $(3,1)$ and $(1, 3)$, the so-called bulge. The subsequent rotations simply chase the bulge out of the lower right corner of the matrix. The rotation in the $(2,3)$ plane sets the elements $(3,1)$ and $(1,3)$ back to zero, but it generates two new nonzero elements at positions $(4,2)$ and $(2,4)$; the rotation in the $(3,4)$ plane sets the elements $(4,2)$ and $(2,4)$ back to zero, but it generates two new nonzero elements at positions $(5,3)$ and $(3,5)$, etc.
Implicit $Q$ Theorem. The effect of this procedure is the following. At the end of the first sweep, the resulting matrix $T_1$ is equal to the the matrix that would have been obtained by factorizing $T-\mu I=QR$ and computing $T_1=RQ+\mu I$.
Since the convergence of the function myTridEigQR() is quadratic (or even cubic),
an eigenvalue is isolated after just a few steps, which requires
$O(n)$ operations. This means that $O(n^2)$ operations are needed to
compute all eigenvalues.
If the eigenvector matrix $Q$ is desired, the plane rotations need to be
accumulated similarly to the accumulation of $X$ in the function myTridiagG().
This accumulation requires $O(n^3)$ operations. Another, faster, algorithm to
first compute only $\Lambda$ and then compute $Q$ using inverse iterations.
Inverse iterations on a tridiagonal matrix are implemented in the LAPACK routine
DSTEIN.
Error bounds. Let $U\Lambda U^T$ and $\tilde U \tilde \Lambda \tilde U^T$ be the exact and the computed EVDs of $A$, respectively, such that the diagonals of $\Lambda$ and $\tilde \Lambda$ are in the same order. Numerical methods generally compute the EVD with the errors bounded by $$ |\lambda_i-\tilde \lambda_i|\leq \phi \epsilon\|A\|_2, \qquad \|u_i-\tilde u_i\|_2\leq \psi\epsilon \frac{\|A\|_2} {\min_{j\neq i} |\lambda_i-\tilde \lambda_j|}, $$ where $\epsilon$ is machine precision and $\phi$ and $\psi$ are slowly growing polynomial functions of $n$ which depend upon the algorithm used (typically $O(n)$ or $O(n^2)$). Such bounds are obtained by combining perturbation bounds with the floating-point error analysis of the respective algorithms.
The eigenvalue decomposition $T=Q\Lambda Q^T$ computed by myTridEigQR() satisfies the
error bounds from fact 7. with $A$ replaced by $T$ and $U$ replaced by $Q$. The deflation criterion
implies $|T_{i,i+1}|\leq \epsilon \|T\|_F$, which is within these
bounds.
The EVD computed by function mySymEigQR() satisfies the error bounds given in Fact 7.
However, the algorithm tends to perform better on matrices, which are
graded downwards, that is, on matrices that exhibit systematic decrease
in the size of the matrix elements as we move along the diagonal.
For such matrices the tiny eigenvalues can usually be computed with higher
relative accuracy (although counterexamples can be easily constructed).
If the tiny eigenvalues are of interest, it should be checked whether
there exists a symmetric permutation that moves larger elements to the
upper left corner, thus converting the given matrix to the
one that is graded downwards.
The function myTridEigQR() is implemented in the
LAPACK subroutine DSTEQR.
This routine can compute just the
eigenvalues, or both eigenvalues and eigenvectors.
The function mySymEigQR() is Algorithm 5 is implemented in the functions eig(), eigvals() and eigvecs(),
and in the LAPACK routine DSYEV.
To compute only eigenvalues, DSYEV calls DSYTRD and DSTEQR
without the eigenvector option. To compute both eigenvalues and eigenvectors,
DSYEV calls DSYTRD, DORGTR, and DSTEQR with the eigenvector
option.
In [44]:
function myTridEigQR(A1::SymTridiagonal)
A=deepcopy(A1)
n=length(A.dv)
T=Float64
λ=Vector{T}(undef,n)
B=Matrix{T}
if n==1
return map(T,A.dv)
end
if n==2
τ=(A.dv[end-1]-A.dv[end])/2
μ=A.dv[end]-A.ev[end]^2/(τ+sign(τ)*sqrt(τ^2+A.ev[end]^2))
# Only rotation
B=A[1:2,1:2]
G,r=givens(B-μ*I,1,2,1)
B=(G*B)*G'
return diag(B)[1:2]
end
steps=1
k=0
while k==0 && steps<=10
# Shift
τ=(A.dv[end-1]-A.dv[end])/2
μ=A.dv[end]-A.ev[end]^2/(τ+sign(τ)*sqrt(τ^2+A.ev[end]^2))
# First rotation
B=A[1:3,1:3]
G,r=givens(B-μ*I,1,2,1)
B=(G*B)*G'
A.dv[1:2]=diag(B)[1:2]
A.ev[1:2]=diag(B,-1)
bulge=B[3,1]
# Bulge chasing
for i = 2 : n-2
B=A[i-1:i+2,i-1:i+2]
B[3,1]=bulge
B[1,3]=bulge
G,r=givens(B,2,3,1)
B=(G*B)*G'
A.dv[i:i+1]=diag(B)[2:3]
A.ev[i-1:i+1]=diag(B,-1)
bulge=B[4,2]
end
# Last rotation
B=A[n-2:n,n-2:n]
B[3,1]=bulge
B[1,3]=bulge
G,r=givens(B,2,3,1)
B=(G*B)*G'
A.dv[n-1:n]=diag(B)[2:3]
A.ev[n-2:n-1]=diag(B,-1)
steps+=1
# Deflation criterion
k=findfirst(abs.(A.ev) .< sqrt.(abs.(A.dv[1:n-1].*A.dv[2:n]))*eps(T))
k=k==nothing ? 0 : k
# display(A)
end
λ[1:k]=myTridEigQR(SymTridiagonal(A.dv[1:k],A.ev[1:k-1]))
λ[k+1:n]=myTridEigQR(SymTridiagonal(A.dv[k+1:n],A.ev[k+1:n-1]))
λ
end
Out[44]:
In [45]:
λ=eigvals(T)
Out[45]:
In [46]:
T
Out[46]:
In [47]:
λ₁=myTridEigQR(T)
Out[47]:
In [48]:
(sort(λ)-sort(λ₁))./sort(λ)
Out[48]:
Once the eigenvalues are computed, the eigeenvectors can be efficiently computed with inverse iterations. Inverse iterations for tridiagonal matrices are implemented in the LAPACK routine DSTEIN.
In [49]:
?LAPACK.stein!
Out[49]:
In [50]:
U=LAPACK.stein!(T.dv,T.ev,sort(λ₁))
Out[50]:
In [51]:
# Orthogonality
norm(U'*U-I)
Out[51]:
In [52]:
# Residual
norm(T*U-U*Diagonal(sort(λ₁)))
Out[52]:
In [53]:
# Some timings - n=100, 200, 400
n=400
Tbig=SymTridiagonal(rand(n),rand(n-1))
@time myTridEigQR(Tbig);
@time λbig=eigvals(Tbig);
@time LAPACK.stein!(Tbig.dv,Tbig.ev,λbig);
In [54]:
n=2000
Tbig=SymTridiagonal(rand(n),rand(n-1))
@time λbig=eigvals(Tbig);
@time U=LAPACK.stein!(Tbig.dv,Tbig.ev,λbig);
@time eigen(Tbig);
Alternatively, the rotations in myTridEigQR() can be accumulated to compute the eigenvectors. This is not optimal, but is instructive. We keep the name of the function, using Julia's multiple dispatch feature.
In [55]:
function myTridEigQR(A1::SymTridiagonal,U::Matrix)
# U is either the identity matrix or the output from myTridiagX()
A=deepcopy(A1)
n=length(A.dv)
T=Float64
λ=Vector{T}(undef,n)
B=Matrix{T}
if n==1
return map(T,A.dv), U
end
if n==2
τ=(A.dv[end-1]-A.dv[end])/2
μ=A.dv[end]-A.ev[end]^2/(τ+sign(τ)*sqrt(τ^2+A.ev[end]^2))
# Only rotation
B=A[1:2,1:2]
G,r=givens(B-μ*I,1,2,1)
B=(G*B)*G'
U*=G'
return diag(B)[1:2], U
end
steps=1
k=0
while k==0 && steps<=10
# Shift
τ=(A.dv[end-1]-A.dv[end])/2
μ=A.dv[end]-A.ev[end]^2/(τ+sign(τ)*sqrt(τ^2+A.ev[end]^2))
# First rotation
B=A[1:3,1:3]
G,r=givens(B-μ*I,1,2,1)
B=(G*B)*G'
U[:,1:3]*=G'
A.dv[1:2]=diag(B)[1:2]
A.ev[1:2]=diag(B,-1)
bulge=B[3,1]
# Bulge chasing
for i = 2 : n-2
B=A[i-1:i+2,i-1:i+2]
B[3,1]=bulge
B[1,3]=bulge
G,r=givens(B,2,3,1)
B=(G*B)*G'
U[:,i-1:i+2]=U[:,i-1:i+2]*G'
A.dv[i:i+1]=diag(B)[2:3]
A.ev[i-1:i+1]=diag(B,-1)
bulge=B[4,2]
end
# Last rotation
B=A[n-2:n,n-2:n]
B[3,1]=bulge
B[1,3]=bulge
G,r=givens(B,2,3,1)
B=(G*B)*G'
U[:,n-2:n]*=G'
A.dv[n-1:n]=diag(B)[2:3]
A.ev[n-2:n-1]=diag(B,-1)
steps+=1
# Deflation criterion
k=findfirst(abs.(A.ev) .< sqrt.(abs.(A.dv[1:n-1].*A.dv[2:n]))*eps())
k=k==nothing ? 0 : k
end
λ[1:k], U[:,1:k]=myTridEigQR(SymTridiagonal(A.dv[1:k],A.ev[1:k-1]),U[:,1:k])
λ[k+1:n], U[:,k+1:n]=myTridEigQR(SymTridiagonal(A.dv[k+1:n],A.ev[k+1:n-1]),U[:,k+1:n])
λ, U
end
Out[55]:
In [56]:
λ,U=myTridEigQR(T,Matrix{Float64}(I,size(T)))
λ
Out[56]:
In [57]:
U
Out[57]:
In [58]:
# Orthogonality
norm(U'*U-I)
Out[58]:
In [59]:
# Residual
norm(T*U-U*Diagonal(λ))
Out[59]:
In [60]:
function mySymEigQR(A::Matrix)
T,H=myTridiag(A)
X=myTridiagX(H)
myTridEigQR(T,X)
end
Out[60]:
In [61]:
A
Out[61]:
In [62]:
λ,U=mySymEigQR(map(Float64,A))
λ
Out[62]:
In [63]:
U
Out[63]:
In [64]:
# Orthogonality
norm(U'*U-I)
Out[64]:
In [65]:
# Residual
norm(A*U-U*Diagonal(λ))
Out[65]:
The $QR$ iterations for unsymmetric matrices are implemented as follows:
The algorithm requires of $O(n^3)$ operations. For more details, see D. Watkins, Unsymmetric Matrix Eigenvalue Techniques and the references therein.
In [66]:
A=rand(-9:9,5,5)
Out[66]:
In [67]:
λ,X=eigen(A)
λ
Out[67]:
In [68]:
X
Out[68]:
In [69]:
# Residual
norm(A*X-X*Diagonal(λ))
Out[69]:
In [70]:
?hessenberg
Out[70]:
In [71]:
B=hessenberg(A)
Out[71]:
In [72]:
B.H
Out[72]:
In [73]:
B.Q
Out[73]:
In [74]:
B.Q'*A*B.Q
Out[74]:
In [77]:
eigvals(Matrix(B.H))
Out[77]:
In [ ]: